Examen final Modules 4 et 5
DUBii 2021
Vichita Ok, d’après Olivier Rué et Valentin Loux
14 April, 2021
Consignes
Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Introduction
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
Analyses
Organisation de votre espace de travail
# répertoire de travail
pwd
# /shared/projects/vok2/m4m5-evaluation
# création de dossiers dans ce répertoire de travail
mkdir analysis
cd analysis
mkdir FASTQ QC CLEANING MAPPING
Téléchargement des données brutes
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
# changement de répertoire
cd..
# pour chercher le nom exact du module à charger (et voir aussi les versions disponibles)
module avail sra-tools
# charger le module
module load sra-tools/2.10.3
# récupération des fichiers fastq du run SRR10390685 dans le dossier FASTQ/
srun --cpus-per-task 6 fasterq-dump --split-files -p SRR10390685 --outdir analysis/FASTQ
# observation générale des fichiers téléchargés
ls -shl analysis/FASTQ/
# compression des fichiers
srun gzip analysis/FASTQ/*.fastq
# taille et format des fichiers compressés
ls -shl analysis/FASTQ/
Combien de reads sont présents dans les fichiers R1 et R2 ?
# nb de reads dans le fichiers R1 (= nb lignes divisé par 4)
srun zcat analysis/FASTQ/SRR10390685_1.fastq.gz | echo $((`wc -l`/4))
# nb de reads dans le fichiers R2
srun zcat analysis/FASTQ/SRR10390685_1.fastq.gz | echo $((`wc -l`/4))Les fichiers FASTQ contiennent 7066055 reads.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
# création d'un dossier REFERENCE
mkdir analysis/REFERENCE
# téléchargement du génome de référence au moyen du lien
cd analysis/REFERENCE
srun wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
# taille et format du fichier téléchargé
ls -lsh
Quelle est la taille de ce génome ?
# format .fasta, donc enlever en premier lieu la première ligne commençant par ">"
# /!\ output de wc: nb lignes | nb mots | nb total caractères
# /!\ enlever les retour à la ligne (= nb de lignes) i.e. soustraire colonne3 - colonne1
srun zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | grep -v ">" | wc | awk '{print $3-$1}'La taille de ce génome est de 4215606 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
srun wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
# taille et format du fichier téléchargé
ls -lshCombien de gènes sont connus pour ce génome ?
# les gènes sont définis par le feature "gene" (colonne3) dans le fichier tabulé
srun zcat GCF_000009045.1_ASM904v1_genomic.gff.gz | grep -c $'\tgene\t'4448 gènes sont recensés dans le fichier d’annotation.
Contrôle qualité
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
# changement de répertoire
cd ../..
# chargement du module fastqc
module avail fastqc
module load fastqc/0.11.9
# fastqc, résultats dans dossier QC
srun --cpus-per-task 8 fastqc analysis/FASTQ/SRR10390685_1.fastq.gz --outdir analysis/QC/ --threads 8
srun --cpus-per-task 8 fastqc analysis/FASTQ/SRR10390685_2.fastq.gz --outdir analysis/QC/ --threads 8
# pour affichage, copier en local, puis visualiser avec navigateur
scp vok2@@core.cluster.france-bioinformatique.fr:/shared/ifbstor1/projects/dubii2021/vok2/m4m5-evaluation/analysis/QC/*.html .
firefox SRR10390685_1_fastqc.html
firefox SRR10390685_2_fastqc.htmlOn peut ainsi vérifier que l’on retrouve les mêmes nombres de reads calculés précédemment.
La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
- Oui
- Non
car - le Phred score est assez élevé (la plupart supérieur à 30) comme le montrent les graphes “Per base sequence quality” et “Per sequence quality scores” - %A=%T et %G=%C et ce pourcentage de 43% en GC est concordant avec l’espèce Bacillus subtilis comme le montre le graphe “Per base sequence content” - aucune position avec “N” comme le montre le graphe “Per base N content”
# chargement du module multiqc
module avail multiqc
module load multiqc/1.9
# multiqc
srun multiqc analysis/QC/
Lien vers le rapport FASTQC R1
Lien vers le rapport FASTQC R2
Lien vers le rapport MulitQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
- Oui
- Non
car tous les reads n’ont pas la même longueur
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
# chargement module seqkit
module load seqkit/0.14.0
# stats pour obtenir le nb de reads et la moyenne
srun seqkit stats --threads 1 analysis/FASTQ/*.fastq.gz
#file format type num_seqs sum_len min_len avg_len max_len
#analysis/FASTQ/SRR10390685_1.fastq.gz FASTQ DNA 7,066,055 1,056,334,498 35 149.5 151
#analysis/FASTQ/SRR10390685_2.fastq.gz FASTQ DNA 7,066,055 1,062,807,718 130 150.4 151
# donc la profondeur de séquençage est égale à :
echo "(7066055*149.5 + 7066055*150.4)/4215606" | bc
La profondeur de séquençage est de : 502 X.
Nettoyage des reads
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
# chargement module fastp
module load fastp/0.20.0
# fastp avec les critères décrits ci-infra
srun --cpus-per-task 8 fastp --in1 analysis/FASTQ/SRR10390685_1.fastq.gz --in2 analysis/FASTQ/SRR10390685_2.fastq.gz --out1 analysis/CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 analysis/CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html analysis/CLEANING/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json /dev/nullLes paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
| length_required | 100 | Taille minimale des reads |
| cut_mean_quality | 30 | Seuil de qualité pour la sélection des reads |
| cut_window_size | 8 | Taille de la fênetre glissante |
| cut_tail | dans le sens 3’-5’ |
# reads restants
srun seqkit stats --threads 1 analysis/CLEANING/*.fastq.gz
#file format type num_seqs sum_len min_len avg_len max_len
#analysis/CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz FASTQ DNA 6,777,048 996,891,051 100 147.1 151
#analysis/CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz FASTQ DNA 6,777,048 990,442,597 100 146.1 151
# pourcentage perte = (initial-final)*100/initial
echo "scale=3;((7066055-6777048)/7006055)*100" | bc
#4.100
Ces paramètres ont permis de conserver 6 777 048 reads pairés, soit une perte de 4.1% des reads bruts.
Alignement des reads sur le génome de référence
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
# chargement module bwa
module load bwa/0.7.17
# alignement
srun --cpus-per-task=16 bwa mem analysis/REFERENCE/GCF_000009045.1_ASM904v1_genomic.fna.gz analysis/CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz analysis/CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 16 > analysis/MAPPING/SRR10390685_on_ASM904v1.sam
# Chargement du module samtools
module load samtools/1.10
# Conversion .sam en .bam
srun --cpus-per-task=8 samtools view --threads 8 analysis/MAPPING/SRR10390685_on_ASM904v1.sam -b > analysis/MAPPING/SRR10390685_on_ASM904v1.bam
# Tri du fichier .bam
srun samtools sort analysis/MAPPING/SRR10390685_on_ASM904v1.bam -o analysis/MAPPING/SRR10390685_on_ASM904v1.sorted.bam
# Indexation du fichier trié .sorted.bam
srun samtools index analysis/MAPPING/SRR10390685_on_ASM904v1.sorted.bam
Combien de reads ne sont pas mappés ?
# Les reads non-mappés correspondent aux flag 4 du fichier .sam.
samtools view -S -f 4 analysis/MAPPING/SRR10390685_on_ASM904v1.sam | wc -l
#744540744540 reads ne sont pas mappés.
Croisement de données
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
# chargement module bedtools
module load bedtools/2.29.2
# Sélection de la ligne du fichier .gff correspondant au gène d'intérêt trmNF (Rq: noms des gènes dans la colonne 3 du fichier .gff)
srun zgrep trmNF analysis/REFERENCE/GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '$3=="gene"' > analysis/REFERENCE/trmNF_gene.gff3
# Recherche des reads qui chevauchent avec au moins 50% (option -f 0.5) de leur longueur le gène trmNF == rechercher le nb de reads mappés sur au moins 50% de leur longueur
srun bedtools intersect -a analysis/MAPPING/SRR10390685_on_ASM904v1.sorted.bam -b analysis/REFERENCE/trmNF_gene.gff3 -f 0.5 > analysis/MAPPING/SRR10390685_on_trmNF_gene.bam
# Tri
srun samtools sort analysis/MAPPING/SRR10390685_on_trmNF_gene.bam -o analysis/MAPPING/SRR10390685_on_trmNF_gene.sorted.bam
# Indexation
srun samtools index analysis/MAPPING/SRR10390685_on_trmNF_gene.sorted.bam
# Regarder le nombre de reads mappés (3ème colonne)
samtools idxstat analysis/MAPPING/SRR10390685_on_trmNF_gene.sorted.bam
#NC_000964.3 4215606 2801 0
#* 0 0 0
2801 reads chevauchent le gène d’intérêt.
Visualisation
Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
# Création du fasta index .fai de notre génome de référence
gunzip analysis/REFERENCE/GCF_000009045.1_ASM904v1_genomic.fna.gz
samtools faidx analysis/REFERENCE/GCF_000009045.1_ASM904v1_genomic.fna
# copier les fichiers du génome (fasta et fasta index) et de l'alignement (bam et bai) dans un répertoire local
capture d’écran
References
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France
